 function results=nwest_balanced_panel(yyy,x,nlag,T,N,lag,~,gls,rep,time_fe,fe) 
% PURPOSE: computes Newey-West adjusted heteroscedastic-serial 
%          consistent Least-squares Regression 
%--------------------------------------------------- 
% USAGE: results = nwest(yyy,x,nlag) 
% where: yyy = dependent variable vector (nobs x 1) 
%        x = independent variables matrix (nobs x nvar) 
%     nlag = lag length to use 
%--------------------------------------------------- 
% RETURNS: a structure 
%        results.meth  = 'newlyw' 
%        results.beta  = bhat 
%        results.tstat = t-stats 
%        results.yhat  = yhat 
%        results.resid = residuals 
%        results.sige  = e'*e/(n-k) 
%        results.rsqr  = rsquared 
%        results.rbar  = rbar-squared 
%        results.dw    = Durbin-Watson Statistic 
%        results.nobs  = nobs 
%        results.nvar  = nvars 
%        results.yyy     = yyy data vector 
% -------------------------------------------------- 
% SEE ALSO: nwest_d, prt(results), plt(results) 
%--------------------------------------------------- 
% References:  Gallant, R. (1987), 
%  "Nonlinear Statistical Models," pp.137-139. 
%--------------------------------------------------- 
 
% written by: 
% James P. LeSage, Dept of Economics 
% University of Toledo 
% 2801 W. Bancroft St, 
% Toledo, OH 43606 
% % jlesage@spatial-econometrics.com 
 
 
 
%if (nargin ~= 3); error('Wrong # of arguments to nwest'); end; 
% clear;
% N=12;T=72;
%  yyy=randn(T*N,1);
% x=randn(T*N,2);
% x(101,:)=nan;yyy(101,:)=nan;
uz=reshape(yyy,T,N);%figuring out different cross section data avaialble for each period;
%N1=sum(~isnan(uz(:,:)'))';%number of avaialble cross section for each period;
[nobs nvar] = size(x); 
 nobs=nobs/N;
results.meth    = 'nwest'; 
results.yyy       = yyy; 
results.nobs    = nobs; 
results.nvar    = nvar; 
 %xx=x;yy=yyy;
%  x(isnan(x)) = 0;yyy(isnan(yyy)) = 0;%replace nan with zeroes
%  xpxi = pinv(x'*x);
%    xpxi = eye(nvar)/(x'*x); 
if time_fe==0
        xpxi = (x'*x)\eye(nvar);
%          xpxi=full(xpxi);
if fe==0
xpxi=sparse(xpxi);
end
else
xpxi = pinv(x'*x);
end    

%         if any(any(isnan(xpxi)))||any(any(xpxi==inf))||any(any(xpxi==-inf))
% %         xpxi = (x'*x+eye(nvar)*10^(-10))\eye(nvar);
% x=full(x);
% xpxi = pinv(x'*x);
%         end
% xpxi=sparse(xpxi);
%            xpxi = pinv(x'*x); %xpxi

% xpxi= ignoreNaN(x'*x,@inv,size(x'*x,2));

results.beta    = xpxi*(x'*yyy); 
results.yhat    = x*results.beta; 
results.resid   = yyy - results.yhat;
% results.resid(isnan(results.resid))=0;

sigu = results.resid'*results.resid; 
results.sigu=sigu;
results.sige    = sigu/(nobs-nvar); 
 if gls==1
     results.sige_gls=reshape(results.resid,length(results.resid)/N,N)'*reshape(results.resid,length(results.resid)/N,N);
results.sige    = sigu; 
 end
uuu=results.resid';% results.goro_fev=sum(results.beta(end-nlag:end,:).^2)*
% perform Newey-West correction 
% emat = []; 
% parfor i=1:nvar; 
% emat = [emat 
%         uuu]; 
% end; 
% emat=zeros(size(uuu,1)*nvar,size(uuu,2));
% BEGINNING OF G_FUNCTION
%     emat=sparse(repmat(uuu,nvar,1));  
%     hhat=sparse(emat.*x'); clear emat x
%     G=zeros(nvar,nvar); w=zeros(2*nlag+1,1); 
%     a=0; 
% %     parfor i=1:nvar
% % u=reshape(hhat(i,:),T,N);
% % % u1=u';
% % uu(i,:)=sum(u,2);
% %     end
% hhat=NDSparse(hhat);
% u=reshape(hhat,T,N,size(hhat,1));
% uu=sum(u,2);clear u hhat
% uu=reshape(uu,nvar,numel(uu)/nvar);
%     while a~=nlag+1; 
%         ga=zeros(nvar,nvar); 
%         w(nlag+1+a,1)=(nlag+1-a)/(nlag+1); 
%         za=uu(:,(a+1):nobs)*uu(:,1:nobs-a)'; 
%           if a==0; 
%            ga=ga+za; 
%           else 
%            ga=ga+za+za'; 
%           end; clear uu
%         G=G+w(nlag+1+a,1)*ga; 
%         a=a+1; 
%     end; % end of while 
%      END OF G_FUNCTION
G=G_FUNCTION(uuu,x,nvar,nlag,T,N,nobs);
%          V1=xpxi*G*xpxi;
% WHAT COMES BELOW (3 LINES) IS BAYESIAN ESTIMATIO REPLACING ABOVE ROW
V=trace(xpxi*G);
%  s=results.resid;
% save('s','s');
%  length(results.resid(results.resid~=0))
% if isnan(sigu)
%     sigu=0.0000001;
% end
% if isnan(V)
%    V=0.0000001;
% end
%  res=results.resid;
% save('xpxi1','xpxi','x','sigu','V','res','yyy');
SIG = iwishrnd(sigu+V,length(results.resid(results.resid~=0)));
          results.beta=mvnrnd(results.beta,1*SIG*xpxi)';



        nwerr= sqrt(diag(V)); 
%         nlag=9;
% if no_interaction==0
% if gls==0
%  results.sig_diff=V1(end-lag+1,end-lag+1)+V1(end-2*lag+1,end-2*lag+1)-2*V1(end-lag+1,end-2*lag+1);
% else
%     results.sig_diff=V(2,2)+V(1,1)-2*V(2,1);
% end  
% % else
%   results.sig_diff_flows=V1(2,2)+V1(1,1)-2*V1(2,1);
% % end   
results.tstat = results.beta./nwerr; % Newey-West t-statistics 
ym = yyy - ones(nobs*N,1)*mean(yyy); 
rsqr1 = sigu; 
rsqr2 = ym'*ym; 
results.rsqr = 1.0 - rsqr1/rsqr2; % r-squared 
%rsqr1 = rsqr1/(nobs-nvar); 
rsqr2 = rsqr2/(nobs-1.0); 
results.rbar = 1 - (rsqr1/rsqr2); % rbar-squared 
ediff = results.resid(2:nobs) - results.resid(1:nobs-1); 
results.dw = diag((ediff'*ediff)./(sigu))'; % durbin-watson 
results.se=nwerr;

% results.se=results.se-imag(results.se);

 if gls==1
     results.se_gls=pinv(x'*x)*sigu;
%              results.se_gls=results.se_gls-imag(results.se_gls);
 end